Exact solutions of the Gross-Pitaevskii equation for stable vortex modes 
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We construct exact solutions of the Gross-Pitaevskii equation for solitary vortices, and approx- 
imate ones for fundamental solitons, in 2D models of Bose-Einstein condensates with a spatially 
modulated nonlinearity of either sign and a harmonic trapping potential. The number of vortex- 
soliton (VS) modes is determined by the discrete energy spectrum of a related linear Schrodinger 
equation. The VS families in the system with the attractive and repulsive nonlinearity are mutually 
complementary. Stable VSs with vorticity S > 2 and those corresponding to higher-order radial 
states are reported for the first time, in the case of the attraction and repulsion, respectively. 
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The realization of Bose-Einstein condensates (BECs) 
in dilute quantum gases has drawn a great deal of in- 
terest to the dynamics of nonlinear excitations in matter 
waves, such as dark H and bright solitons vortices 
[11,13], supervortices Q, etc. The work in this direction 
was strongly stimulated by many similarities to solitons 
and vortices in optics 6]. In the mean-field approx- 
imation, the BEC dynamics at ultra-low temperatures 
is accurately described by the Gross-Pitaevskii equation 
(GPE), the nonlinearity being determined by the s-wave 
scattering length of interatomic collisions, which can be 
controlled by means of the magnetic 0] or low-loss opti- 
cal 0] Feshbach-resonance technique, making spatiotem- 
poral "management" of the local nonlinearity possible 
through the use of time-dependent and/or non- uniform 
fields p. In particular, a recently developed technique, 
which allows one to "paint" arbitrary one- and two- 
dimensional (ID and 2D) average patterns by a rapidly 
moving laser beam suggests new perspectives for the 
application of the nonlinearity management based on the 
optical Feshbach-resonance. A counterpart of the GPE, 
which is a basic model in nonlinear optics, is the nonlin- 
ear Schrodinger equation (NLSE) [6]. In the latter case, 
the modulation of the local nonlinearity can be imple- 
mented too-in particular, by means of indiffusion of a 
dopant resonantly interacting with the light. 

The search for exact solutions to the GPE /NLSE is an 
essential direction in the studies of nonlinear matter and 
photonic waves. In particular, stable soliton solutions are 
of great significance to experiments and potential appli- 
cations, as they precisely predict conditions which should 
allow the creation of matter- wave and optical solitons, 
including challenging situations which have not yet been 
tackled in the experiment, such as 2D matter-wave soli- 
tons and solitary vortices [lO(. In the ID setting with the 
spatially uniform nonlinearity, exact stable solutions for 
bright and dark solitons have been constructed in special 
cases, when the external potential is linear or quadratic 



[IH, and a specially designed inhomogeneous nonlinear- 
ity may support bound states of an arbitrary number of 



solitons 12[. In the 2D geometry with the uniform non- 



linearity, exact delocalized solutions were constructed in 
the presence of a periodic potential [l3| . In the 2D geom- 
etry, ID bright and dark solitons are unstable, with the 
bright ones suffering the breakup into a chain of collaps- 
ing pulses, and the dark solitons splitting into vortex- 
antivortex pairs [J [13] • Another approach to finding ex- 
act solitonic states in both ID and 2D settings (with the 
uniform nonlinearity), along with the supporting poten- 
tials, was proposed in the form the respective inverse 
problem |15j. As concerns vortices, while delocalized 
ones have been created and extensively investigated in 
self-repulsive BECs and vortex solitons (VSs) were 
created in photorefractive crystals equipped with pho- 
tonic lattices [16| . analytical solutions for localized vor- 
tices have not been reported yet, and experimental cre- 
ation of stable VSs in self-attractive BEC and optical 
media with the fundamental cubic nonlinearity remains 
a challenge to the experiment. 

In this work we demonstrate that exact VS solutions 
can be found in the 2D GPE/NLSE, with a specially de- 
signed (and experimentally realizable) profile of the ra- 
dial modulation of the nonlinearity coefficient, of both 
attractive and repulsive types. In fact, the search for lo- 
calized solutions in 2D models with diverse nonlinearity- 
modulation profiles has recently attracted considerable 
attention [l7(, although no exact solutions have been re- 
ported. We consider only nonlinearity-modulation pat- 
terns which do not change the sign, as zero-crossing 
would make the use of the FR problematic 0, H| • 

The scaled form of the GPE/NLSE is iip t = -V 2 i/; + 
g(r)\i/j\ 2 4'+V(' r )' t Pi where tp is the BEC macroscopic wave 
function, V 2 the 2D Laplacian, and g(r) is the nonlinear- 
ity coefficient which, as well as external potential V(r), 
is a function of radial coordinate r. The same equa- 
tion finds a straightforward implementation as a spatial- 
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FIG. 1: (Color online), (a) Exact vortex solitons without the 
external potential. Inset, the corresponding profiles of the 
attractive-nonlinearity coefficient, (c) Exact vortex solitons of 
different radial quantum number n with S = 1, the respective 
nonlinearity-coemcient profile being depicted by the solid line 
in the inset of (a). Parameters are 03,4 = — p. = —go = 1. (b) 
and (d): The same as (a) and (c) when the trapping potential 
is present, V = O.Olr 2 , with ci,2 = 3. 

domain model of the beam propagation in bulk media, 
with t, ip, and g(r) being the propagation distance, am- 
plitude of the electromagnetic field, and the local Kerr 
coefficient. 

Assuming the stationary wave function as ip{r, 9, t) = 
4>{r) exp(iS6 — ip£), where 9 is the azimuthal angle, S an 
integer vorticity, and p the chemical potential (or prop- 
agation constant in optics) , leads to the equation for the 
real stationary wave function 4>(r): pip = — 0" — r~ l (f)' + 
g{r)(j> 3 + [S 2 r~ 2 +V(r)\ <j>. For S 7^ 0, <j>(r) should vary 
as r' 5 ! at r — > 0, which is replaced by <p'(r = 0) = for 
S = 0. The localization requires <p(r = 00) = 0. 

Defining 0(r) = p(r)U(R(r)), g(r) = go r ~ 2 P~ 6 { r ), 
with R(r) = Jq s~ 1 p~ 2 (s)ds, one can find that p(r) and 
U(R) obey the following equations: 

p" + r- l p' + (p- V(r) - S 2 r- 2 ) p = Er^p' 3 , (I) 
-d 2 U/dR 2 + g U 3 = EU, (2) 

where E and go are constants. The reduction of the 
GPE/NLSE to Eq. ^ helps one to find exact solutions, 
as the latter equation is solvable in terms of elliptic func- 
tions. Then, if a solution to Eq. (TTJ) is known, one can 
construct exact solutions to the underlying GPE/NLSE, 
while the nonlinearity-modulation profile admitting ex- 
act solutions is determined by p(r). Physical solutions 
impose restrictions on p: from expressions for R(r) and 
g(r) it follows that p cannot change its sign; further, it 
must behave as r~ a with a > 1/3 at r - > 0, and diverge 
(p — > 00) at r — > 00, so that the nonlinearity strength is 
bounded and the integration in R(r) converges. 

We begin constructing exact VS solutions for the at- 
tractive nonlinearity (go < 0) when E — 0, so that Eq. 
(TTJ) is solvable. With the harmonic potential V = kr 2 , 
p can be found in terms of the Whittaker's M and W 
functions [18]: p(r) = r" 1 [ Cl M(p/4Vk, \S\/2, Vkr 2 ) + 
C2W(/i/4v / fc, 15*1/2, Vkr 2 )], where the restrictions on p re- 
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FIG. 2: (Color online), (a) The number of numerically found 
vortex-soliton modes versus the chemical potential in the case 
of the repulsive nonlinearity, with the harmonic trap, (b) The 
largest instability growth rate for numerically found vortices 
with S = 2, E = 1 and k = 0.01. 

quire p < po = 2(1 + I^Dn/A: and C1C2 > 0. Without the 
trap (k = 0), p degenerates to p(r) = c^Isi^—pr) + 
CAKs{\J—pr), with p < po — 0, Is and K$ being the 
modified Bessel and Hankel functions, and the constants 
satisfying C3C4 > 0. In both cases, one has p(r) ~ r~l s l 
at r — > 0, hence g(r) ~ r 6 l s l~ 2 and R(r) ~ r 2 l s l at r 0, 
and p(r) —¥ 00 as r — > 00. Thus, the respective nonlin- 
earity is localized and bounded, and R{r) is bounded too. 
To meet boundary conditions 0(0) = </>(oo) = 0, an exact 
solution to Eq. ([2]) is chosen as 

U{R) = (n»/A/=flo) cn (nrjR - K(l/y/2), I/V2) , (3) 

where n = 2, 4, 6, • • • with n/2 being the radial quantum 
number, rj = K{y/2/2)/R{r = 00), and K(l/Vf) is the 
complete elliptic integral of the first kind. 

It follows from Eq. © that U(R) ~ R at R -t 0, which 
implies that the amplitude of the exact VS is pU ~ r' 5 ' 
at r — > 0, as it should be. Thus, for given p, S 1 , and non- 
linearity strength go (and fc, in the presence of the trap), 
one can construct an infinite number of exact VSs with 
n/2 bright rings surrounding the vortex core, as shown 
in Fig. 1. Although these VSs share the same chemi- 
cal potential, their energies increase with the increase of 
even number n. 

Next we consider the repulsive nonlinearity, go > 0, 
and include the harmonic trap to confine the system. In 
this case, the existence of elliptic-function solutions to 
Eq. @ requires E > 0, hence Eq. |T]) is a nonlinear 
equation, which can be solved only in a numerical form. 
To construct the VSs in this case, one requires p ~ r~' s ' 
(for 5 7^ 0) at r — > 0, hence the nonlinear term in Eq. 
(TT]) may be neglected near r = 0. Thus, p is similar to 
the Neumann function, Ys(y/Jlr), at r — > for p > (for 
p < it can be checked that VS solutions do not exist). 
On the other hand, p — > 00 at r — > 00, due to the the 
presence of the harmonic trap. Further, term Er~ 2 p~ 3 
with E > in Eq. (TTJ) guarantees the sign-definiteness of 
p(r). Therefore, taking small ro as an initial point and 
using the Neumann function and its derivative at r = ?*o 
as initial conditions, one can numerically integrate Eq. 



FIG. 3: (Color online), (a) Modulation profiles of the re- 
pulsive nonlinearity for numerically found vortices, which are 
displayed in (b) for S = 1 and in (c) for S = 2, in the presence 
of the harmonic trap. Opened circles, squares, solid circles, 
and triangles denote solutions with n = 2, 4, 6, and 8, respec- 
tively. Parameters are E = 1, fi — 2, k — go = 0.01. 

([T]) to obtain R(r) and <?(r); then VSs can be constructed 
in the numerical form, using the exact solution to Eq. ([2} , 

U(R) = ^2{E- B 2 )/g an (bR, ^E/B 2 - l) , (4) 

subject to a constraint with even numbers n, 

y/W/2 <B = nK(\J E I B 2 - i)/R(ao) < Ve. (5) 

From Eq. flS)) it follows that n < n max = 2R(oo)VE/n. 
This implies that there is a finite number of the VS modes 
(or none, if n max < 2), in contrast to the case of the at- 
tractive nonlinearity, where the infinite set of exact VSs 
was constructed. Figure 2(a) shows the number of the 
numerical VS solutions versus /i, demonstrating that the 
cutoff value of the chemical potential is same as for the 
exact VSs with the attractive nonlinearity, and the num- 
ber of VSs jumps at points p — \if' = 2(2j + |5| — l)Vk, 
with j = 1, 2, 3, • • • , which is precisely the j-th energy 
eigenvalue of the vortex state in the corresponding lin- 
ear Schrodinger equation. Thus, there are j numeri- 
cally exact VSs, associated with expression Q, where 

(S) (S) 

n = 2,4, 6, • • • , 2j, in interval p) < p < p)+i- A simi- 
lar result is true for the ID GPE with the optical-lattice 
potential, where n families of fundamental gap solitons 
exist in the n-th bandgap fl9j ] . Figure 3 displays a char- 
acteristic example of the numerically found VSs, in the 
case when four of them exist. Similarly, it can be shown 
that there are infinitely many VSs, but just the first j 
modes do not exist, when Pj<p< /uf+i for the attrac- 
tive nonlinearity, which demonstrates that the repulsive 
and attractive nonlinearities are mutually complementary 
ones, in this respect. 

Solutions can also be found for 2D fundamental soli- 
tons (S — 0) supported by the following two-tier nonlin- 
earity, with constant g r : 

g(r) = g r (at < r < r ), g Q p~ 6 r~ 2 (at r > r ) . (6) 

For r > ro, exact solutions can be found in the same 
way as above, except that now R(r) = ds/[s (p(s)) 2 }. 
At r < ro, if ro is small in comparison with the spatial 
scale of the external potential (for instance, ro -C y/TJk, 
in the presence of V = kr 2 ), 4>(r) may be approxi- 
mated by a constant, cj> = — V(0)]/g r . Because 4>(r) 



FIG. 4: (Color online). Fundamental solitons supported by 
the repulsive two-tier nonlinearity, see Eq. JB). The con- 
stant values of g r for n — 1,3,5,7 and 9 are respectively 
0.0200,0.0201,0.0224,0.0325, and 0.1093. Parameters are 
p{r ) = 1, p'(r ) =0,r = E = l,k = g = 0.01, fx = 2. 



and 4>'(r) must be continuous at r = ro, one then re- 
quires p'(r ) = and dU(0)/dR = 0, which leads to 
g r = [/i — V(0)] [p(r )U(0)]~ 2 . In this case, the solution 
to Eq. © is given by Eqs. Q and (P, with BR re- 
placed by BR + K(y/E/B 2 - 1), and n = 1, 3, 5, ■ ■ • for 
the repulsive nonlinearity, to make <ft(oo) — 0. Similar 
solutions can be constructed for the attractive nonlinear- 
ity. Examples of the fundamental solutions are presented 
in Fig. 4. 

We employed the linear stability analysis and direct 
simulations to verify the stability of the solitons. For 
the attractive nonlinearity, we have found that all exact 
VSs are subject to the azimuthal modulational instability 
without the external potential. As a result, the VSs break 
up, and eventually collapse. However, the lowest-order 
exact VSs with n — 2 are stable in the presence of the 
harmonic trap, if p is close enough to the cutoff value po, 
while the amplitude of the VS does not become small, see 
Fig. 5. We stress that, in numerous previous works, a 
stability region for VSs trapped in the harmonic potential 
was found solely for S = 1, all vortices with S > 2 being 
conjectured unstable The present results report the 
first example of stable vortices with S > 2 in the trapped 
self- attractive fields. 

In the case of the repulsive nonlinearity combined with 
the harmonic trap, the numerically found VSs can be sta- 
ble for every n at which they exist, within some region 
of values of p, see Fig. 2(b). This finding is remarkable 
in the sense that previous works did not report stable 
trapped vortices in higher-order radial states. Moreover, 
Figs. 2(b) and 6 demonstrates that VSs with the largest 
number of rings may be more stable than its counter- 
parts with fewer rings. An explanation to this feature 
is that local density peaks place themselves in troughs 
of the nonlinearity landscape, thus lowering the system's 
energy. Similar properties are featured by the fundamen- 
tal solitons. Numerical simulations show that those VSs 
with n = 2 which are unstable either split into vortices 
with lower topological charges or exhibit a quasi-stability, 
periodically breaking and recovering the axial symmetry 
[Fig. 6(a)], similar to what was previously observed in 
the case of the attractive nonlinearity [3, [13], while unsta- 
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FIG. 5: (Color online), (a) The largest instability growth rate 
for exact vortex solitons with S = 3 and the attractive non- 
linearity, (b) A stable vortex (solid lines, with circles showing 
its profile at t — 200, when it was initially perturbed by ran- 
dom noise) , and (c) the corresponding nonlinearity coefficient 
when n — 0.7. Here go = —1000, k — 0.01, and ci,2 = 3. 



20 




20 




20 




10 


^^^^ 






10 




>. 
-10 




: 




*Q 
-It) 




-20 




-20 




-20 





20- -10 10 20 -20 -10 10 20 -20 -10 10 20 

XXX 



FIG. 6: (color online), (a) The result of the quasi-stable 
evolution of a numerically found vortex with n — 2, in the 
case of the self-repulsion, at t = 1700. (b, c) Unstable vortices 
with n = 4 and 6 at t — 60 and 70, respectively. Parameters 
are S = 2, E = 1, g = k = 0.01, fx = 1.82. The vortex with 
n = 8 are stable. 

ble VSs with n = 4, 6, • • • ultimately evolve into vortices 
located close to zero-amplitude points. 

In conclusion, we have constructed exact solutions of 
the Gross-Pitaevskii equation for solitary vortices, and 
approximate ones for fundamental solitons in the frame- 
work of the GPE/NLSE with 2D axisymmetric profiles of 
the modulation of the nonlinearity coefficient, and har- 
monic trapping potential. We have demonstrated that 
the attractive/repulsive nonlinearity supports an infi- 
nite/finite number of exact VSs. In particular, stable 
VSs with vorticity S > 2, as well as those correspond- 
ing to higher-order radial states, have been produced 
for the first time. The results suggest a scenario for 
the creation of stable vortex solitons in BEC and optics, 
which have not been as yet observed in experiments. The 
necessary BEC nonlinearity landscape can be built by 
means of the Feshbach-resonance technique. The corre- 
sponding nonuniform magnetic field may be created by a 
micro-fabricated ferromagnetic structure integrated into 
the matter- wave chip [21], or one can use the respective 
pattern of the laser beams. In optics, the same nonlinear- 
ity landscape may created by a nonuniform distribution 
of nonlinearity-enhancing dopants. The method devel- 
oped in this work for finding the exact solutions, which 
is based on reducing the 2D equation to the solvable sys- 
tem of Eqs. (QJ and ([2]), can be applied to other models. 
In particular, a challenging problem is to devise a physi- 
cally relevant model admitting exact 3D solitons. 
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